Heatmaps

Load packages

suppressPackageStartupMessages({
  library(tidyverse)
  library(pheatmap)
  library(here)
  library(ggiraph)
  library(ggrepel)
})

Differential Expression (DTE & DGE) Heatmaps

Load Transcript CPM & DGE data

Load gene CPM & DGE data

Define Subsets

top_tx_df <- tx_df %>%
  mutate(is_sig = FDR_DEG < 0.05, abs_logFC = abs(logFC_DEG)) %>%
  arrange(desc(is_sig), desc(abs_logFC))

top100_tx <- top_tx_df %>% slice_head(n = 100) %>% pull(row_id)
top50_tx <- head(top100_tx, 50)

top_gene_df <- gene_df %>%
  mutate(is_sig = FDR_DEG < 0.05, abs_logFC = abs(logFC)) %>%
  arrange(desc(is_sig), desc(abs_logFC))

top100_genes <- top_gene_df %>% slice_head(n = 100) %>% pull(row_id)
top50_genes <- head(top100_genes, 50)

Heatmap function for transcript CPM

plot_tx_heatmap <- function(tx_ids, title) {
  mat <- cpm_matrix[rownames(cpm_matrix) %in% tx_ids, , drop = FALSE]
  mat <- log2(mat + 1)
  mat <- mat[rowSums(mat) > 0, , drop = FALSE]
  if (nrow(mat) == 0) return(NULL)

  # Use the modified label map with PB accessions and asterisks
  rownames(mat) <- tx_label_map[rownames(mat)]

  annotation_df <- data.frame(
    sample_name = colnames(mat),
    group = ifelse(grepl("WT", colnames(mat)), "WT", "Q157R")
  ) %>% column_to_rownames("sample_name")

  annotation_colors <- list(group = c(Q157R = "#336B87", WT = "#E99787"))

  col_order <- c(grep("Q157R", colnames(mat), value = TRUE),
                 grep("WT", colnames(mat), value = TRUE))
  mat <- mat[, col_order]
  annotation_df <- annotation_df[col_order, , drop = FALSE]

  n_rows <- nrow(mat)
  plot_height <- min(10, 5 + n_rows * 0.15)
  fontsize_row <- if (n_rows <= 30) 10 else if (n_rows <= 60) 8 else 6

  pheatmap(mat, scale = "row", cluster_rows = TRUE, cluster_cols = FALSE,
           annotation_col = annotation_df, annotation_colors = annotation_colors,
           show_colnames = TRUE, show_rownames = TRUE, main = title,
           fontsize_row = fontsize_row, fontsize_col = 10, height = plot_height)
}

Heatmap function for gene CPM

plot_gene_heatmap <- function(gene_ids, title) {
  mat <- gene_cpm[rownames(gene_cpm) %in% gene_ids, , drop = FALSE]
  mat <- log2(mat + 1)
  mat <- mat[rowSums(mat) > 0, , drop = FALSE]
  if (nrow(mat) == 0) return(NULL)

  rownames(mat) <- gene_label_map[rownames(mat)]

  annotation_df <- data.frame(
    sample_name = colnames(mat),
    group = ifelse(grepl("WT", colnames(mat)), "WT", "Q157R")
  ) %>% column_to_rownames("sample_name")

  annotation_colors <- list(group = c(Q157R = "#336B87", WT = "#E99787"))

  col_order <- c(grep("Q157R", colnames(mat), value = TRUE),
                 grep("WT", colnames(mat), value = TRUE))
  mat <- mat[, col_order]
  annotation_df <- annotation_df[col_order, , drop = FALSE]

  n_rows <- nrow(mat)
  plot_height <- min(10, 5 + n_rows * 0.15)
  fontsize_row <- if (n_rows <= 30) 10 else if (n_rows <= 60) 8 else 6

  pheatmap(mat, scale = "row", cluster_rows = TRUE, cluster_cols = FALSE,
           annotation_col = annotation_df, annotation_colors = annotation_colors,
           show_colnames = TRUE, show_rownames = TRUE, main = title,
           fontsize_row = fontsize_row, fontsize_col = 10, height = plot_height)
}

Plot Top Transcript & Gene Heatmaps

plot_tx_heatmap(top50_tx, "Top 50 Differentially Expressed Transcripts")

plot_tx_heatmap(top100_tx, "Top 100 Differentially Expressed Transcripts")

plot_gene_heatmap(top50_genes, "Top 50 Differentially Expressed Genes")

plot_gene_heatmap(top100_genes, "Top 100 Differentially Expressed Genes")

Differential Transcript Usage (DTU) Heatmaps

Load DTU Data

Define DTU Subsets

top_tx_dtu_df <- tx_dtu_df %>%
  arrange(adj.p.value_DTU, desc(abs(lf_DTU)))

top100_tx_dtu <- top_tx_dtu_df %>% slice_head(n = 100) %>% pull(row_id)
top50_tx_dtu <- head(top100_tx_dtu, 50)

top_gene_dtu_df <- gene_dtu_df %>%
  arrange(adj_p.value_DTU, desc(abs(lf_DTU)))

top100_gene_dtu <- top_gene_dtu_df %>% slice_head(n = 100) %>% pull(row_id)
top50_gene_dtu <- head(top100_gene_dtu, 50)

Heatmap function for DTU – note, it is not common to show DTU with heatmaps, so I am throwing in a volcano plot as well

plot_dtu_heatmap <- function(mat, row_label_map, row_ids, title) {
  mat <- mat[rownames(mat) %in% row_ids, , drop = FALSE]
  mat <- mat[rowSums(mat) > 0, , drop = FALSE]
  mat <- mat[apply(mat, 1, function(x) sd(x, na.rm = TRUE) > 0), , drop = FALSE]
  if (nrow(mat) == 0) return(NULL)

  rownames(mat) <- row_label_map[rownames(mat)]

  annotation_df <- data.frame(
    sample_name = colnames(mat),
    group = ifelse(grepl("WT", colnames(mat)), "WT", "Q157R")
  ) %>% column_to_rownames("sample_name")

  annotation_colors <- list(group = c(Q157R = "#336B87", WT = "#E99787"))

  col_order <- c(grep("Q157R", colnames(mat), value = TRUE),
                 grep("WT", colnames(mat), value = TRUE))
  mat <- mat[, col_order]
  annotation_df <- annotation_df[col_order, , drop = FALSE]

  n_rows <- nrow(mat)
  plot_height <- min(10, 5 + n_rows * 0.15)
  fontsize_row <- if (n_rows <= 30) 10 else if (n_rows <= 60) 8 else 6

  pheatmap(mat, scale = "row", cluster_rows = TRUE, cluster_cols = FALSE,
           annotation_col = annotation_df, annotation_colors = annotation_colors,
           show_colnames = TRUE, show_rownames = TRUE, main = title,
           fontsize_row = fontsize_row, fontsize_col = 10, height = plot_height)
}

Plot DTU Transcript Heatmaps

plot_dtu_heatmap(frac_matrix, frac_label_map, top100_tx_dtu, "Top 100 DTU Transcripts")

plot_dtu_heatmap(frac_matrix, frac_label_map, top50_tx_dtu, "Top 50 DTU Transcripts")

Plot DTU Gene Heatmaps

plot_dtu_heatmap(gene_frac_matrix, gene_frac_label_map, top100_gene_dtu, "Top 100 DTU Genes")

plot_dtu_heatmap(gene_frac_matrix, gene_frac_label_map, top50_gene_dtu, "Top 50 DTU Genes")

Interactive Volcano plot

Alternative Splice Event Heatmaps

Load Alternative Splicing Data

Define Alternative Splicing Subsets

top_as_df <- as_df %>%
  arrange(`Q157R-WT_p.value`, desc(abs(`Q157R-WT_dPSI`)))

top100_as <- top_as_df %>% slice_head(n = 100) %>% pull(row_id)
top50_as <- head(top100_as, 50)

Heatmap function for Alternative Splicing

plot_as_heatmap <- function(mat, row_label_map, row_ids, title) {
  # Filter matrix to include selected row_ids
  mat <- mat[rownames(mat) %in% row_ids, , drop = FALSE]
  
  # Keep rows with non-zero sum and some variability
  mat <- mat[rowSums(mat) > 0, , drop = FALSE]
  mat <- mat[apply(mat, 1, function(x) sd(x, na.rm = TRUE) > 0), , drop = FALSE]
  if (nrow(mat) == 0) return(NULL)

  # Apply human-readable row labels (must be unique)
  rownames(mat) <- row_label_map[rownames(mat)]

  # Column annotation: WT vs Q157R
  annotation_df <- data.frame(
    sample_name = colnames(mat),
    group = ifelse(grepl("WT", colnames(mat)), "WT", "Q157R")
  )
  rownames(annotation_df) <- annotation_df$sample_name
  annotation_df$sample_name <- NULL

  # Annotation colors
  annotation_colors <- list(
    group = c(Q157R = "#336B87", WT = "#E99787")
  )

  # Column order: Q157R first, then WT
  col_order <- c(grep("Q157R", colnames(mat), value = TRUE),
                 grep("WT", colnames(mat), value = TRUE))
  mat <- mat[, col_order]
  annotation_df <- annotation_df[col_order, , drop = FALSE]

  # Dynamic plot sizing
  n_rows <- nrow(mat)
  plot_height <- min(10, 5 + n_rows * 0.15)
  fontsize_row <- if (n_rows <= 30) 10 else if (n_rows <= 60) 8 else 6

  # Final heatmap plot
  pheatmap::pheatmap(
    mat,
    scale = "row",
    cluster_rows = TRUE,
    cluster_cols = FALSE,
    annotation_col = annotation_df,
    annotation_colors = annotation_colors,
    show_colnames = TRUE,
    show_rownames = TRUE,
    main = title,
    fontsize_row = fontsize_row,
    fontsize_col = 10,
    height = plot_height
  )
}

Plot Alternative Splicing Heatmaps - these use a typical internal standard to show change

plot_as_heatmap(as_matrix, as_label_map, top100_as, "Top 100 Splicing Events")

plot_as_heatmap(as_matrix, as_label_map, top50_as, "Top 50 Splicing Events")

## Alternative Splicing heatmaps with dPSI

plot_psi_heatmap <- function(as_df, row_ids, title) {
  # Define PSI columns
  psi_cols <- c("A258_Q157R_PSI", "X504_Q157R_PSI", "A309_Q157R_PSI",
                "V335_WT_PSI", "V334_WT_PSI", "A310_WT_PSI")
  
  # Filter to selected row_ids
  filtered_df <- as_df %>%
    filter(row_id %in% row_ids) %>%
    select(row_id, label, all_of(psi_cols), is_sig)
  
  # Create PSI matrix
  psi_matrix <- filtered_df %>%
    select(row_id, all_of(psi_cols)) %>%
    column_to_rownames("row_id") %>%
    as.matrix()
  
  # Keep rows with non-zero sum and some variability
  psi_matrix <- psi_matrix[rowSums(psi_matrix, na.rm = TRUE) > 0, , drop = FALSE]
  psi_matrix <- psi_matrix[apply(psi_matrix, 1, function(x) sd(x, na.rm = TRUE) > 0), , drop = FALSE]
  
  if (nrow(psi_matrix) == 0) {
    warning("No valid PSI data found for the selected row_ids")
    return(NULL)
  }
  
  # Apply human-readable row labels (ensure uniqueness)
  remaining_row_ids <- rownames(psi_matrix)
  remaining_labels <- filtered_df$label[match(remaining_row_ids, filtered_df$row_id)]
  rownames(psi_matrix) <- make.unique(remaining_labels)
  
  # Clean up column names
  colnames(psi_matrix) <- gsub("_PSI", "", colnames(psi_matrix))
  
  # Column annotation: sample type
  col_annotation_df <- data.frame(
    group = ifelse(grepl("WT", colnames(psi_matrix)), "WT", "Q157R"),
    row.names = colnames(psi_matrix)
  )
  
  # Annotation colors
  annotation_colors <- list(
    group = c(Q157R = "#336B87", WT = "#E99787")
  )
  
  # Column order: Q157R first, then WT
  col_order <- c(grep("Q157R", colnames(psi_matrix), value = TRUE),
                 grep("WT", colnames(psi_matrix), value = TRUE))
  psi_matrix <- psi_matrix[, col_order]
  col_annotation_df <- col_annotation_df[col_order, , drop = FALSE]
  
  # Dynamic plot sizing
  n_rows <- nrow(psi_matrix)
  plot_height <- min(12, 6 + n_rows * 0.2)
  fontsize_row <- if (n_rows <= 30) 10 else if (n_rows <= 60) 8 else 6
  
  # Create heatmap
  pheatmap::pheatmap(
    psi_matrix,
    scale = "row",  # Scale by row to compare patterns across samples
    cluster_rows = TRUE,
    cluster_cols = FALSE,
    annotation_col = col_annotation_df,
    annotation_colors = annotation_colors,
    show_colnames = TRUE,
    show_rownames = TRUE,
    main = title,
    fontsize_row = fontsize_row,
    fontsize_col = 10,
    height = plot_height
  )
}

Plot

plot_psi_heatmap(as_df, top50_as, "Top 50 Alternative Splicing Events")

plot_psi_heatmap(as_df, top100_as, "Top 100 Alternative Splicing Events")